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Abstract 

The task of camera calibration is to estimate the intrinsic and extrinsic parameters of a camera 
model. Though there are some restricted techniques to infer the 3-D information about the scene 
from uncalibrated cameras, effective camera calibration procedures will open up the possibility of 
using a wide range of existing algorithms for 3-D reconstruction and recognition. The applications 
of camera calibration include vision-based metrology, robust visual platooning and visual docking of 
mobile robots where the depth information is important. 

Key Words: Camera calibration, projective model, radial distortion, radial undistortion, flexible 
setup, nonlinear optimization. 
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Table 1: List of Variables 



Variable 



Explanation 



XumlxT of fcatuix^ points in each ima,t!,(> 



N 



Number of images taken for calibration 



T 



Focal length 



Skewness angle between two image axes 



3-D point in world reference frame 



3-D point in camera reference frame 



2-D point in camera frame with z = f 



3-D point in world reference frame with Z'^ = 



The corresponding projected point in image plane of Mj 



a, (3, J, uo,vo 



5 intrinsic parameters 



Effective sizes of the pixel in the horizontal and vertical directions 



(/a: fy 



fx = f/Sx, fy = f/i 



k = (fcl, k2 



Distortion coefficients 
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Real observed distorted image points 



Ideal projected undistorted image points 
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Scaling factor 
Objective function 



J 



A = 



a 7 Mo 
1 



Camera intrinsic matrix 



B = A-^ A-i 



Absolute conic 



R 



.'-) X .'-) rotation matrix 



3x1 translation vector 



[RtT 



Camera extrinsic matrix 



H 



Homography matrix 



Matrix used to estimate homography matrix H 



Matrix stacking constraints to estimate intrinsic parameters 
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1 Introduction 



Depending on what kind of calibration object being used, there are mainly two categories of calibration 
methods: photogrammetric calibration and self-calibration. Photogrammetric calibration refers to those 
methods that observe a calibration object whose geometry in 3-D space is known with a very good 
precision. Self-calibration does not need a 3-D calibration object. Three images of a coplanar object 
taken by the same camera with fixed intrinsic camera parameters are sufficient to estimate both intrinsic 
and extrinsic parameters. The obvious advantage of the self-calibration method is that it is easy to 
set up and the disadvantage is that it is usually considered unreliable. However, the author of [3] 
shows that by preceding the algorithm with a very simple normalization (translation and rotation 
of the coordinates of the matched points), results arc obtained comparable with the best iterative 
algorithms. A four step calibration procedure is proposed in [4]. The four steps are: linear parameter 
estimation, nonlinear optimization, correction using circle/ellipse, and image correction. But for a 
simple start, linear parameter estimation and nonlinear optimization are enough. In [5], a plane-based 
calibration method is described where the calibration is performed by first determining the absolute 
conic B = A^-^A"-*^, where A is a matrix formed by a camera's 5 intrinsic parameters (See Section 4.2). 
In [5], the parameter 7 (a parameter describing the skewness of the two image axes) is assumed to be 
zero and the authors observe that only the relative orientations of planes and camera is of importance for 
singularities and planes that are parallel to each other provide exactly the same information. Recently, 
Intel distributes its "Camera Calibration Toolbox for Matlab" freely available online [6]. The Intel 
camera calibration toolbox first finds the feature locations of the input images, which are captTirod by 
the camera to be calibrated using a checkerboard calibration object. Then, it calculates the camera's 
intrinsic parameters. However, when we used the images captured by our desktop camera as the input 
images, the detected feature locations contain great errors. We decided not to use Intel's method since 
its flexibility and accuracy are poor. Therefore, in this report, our work is mainly based on the self- 
calibration algorithm originally developed by Microsoft Research Group [1, 2], which has been commonly 
regarded as a great contribution to the camera calibration. The key feature of Microsoft's calibration 
method is that the absolute conic B is used to estimate the intrinsic parameters and the parameter 7 is 
considered. The proposed technique in [1, 2] only requires the camera to observe a planar pattern at a 
few (at least 3, if both the intrinsic and the extrinsic parameters are to be estimated uniquely) different 
orientations. Either the camera or the calibration object can be moved by hand as long as they cause 
no singularity problem and the motion of the calibration object or camera itself need not to be known 
in advance. 

By "flexibility", we mean that the calibration object is coplanar and easy to setup while by "robust- 
ness" , it implies that the extracted feature locations are accurate and the possible singularities due to 
improperly input images can be detected and avoided. 

The main contributions in this report are briefly summarized as follows: 

(1) A complete code platform implementing Microsoft's camera calibration algorithm has been built. 

(Microsoft did not release the feature location prc-proccssor [1, 2]); 

(2) A technical error in Microsoft's camera calibration equations has been corrected (Equation (26)); 

(3) A new method to effectively find the feature locations of the calibration object has been used 
in the code. More specifically, a scan line approximation algorithm is proposed to accurately 
determine the partitions of a given set of points; 

(4) A numerical indicator is used to indicate the possible singularities among input images to enhance 
the robustness in camera calibration (under development and to be included); 

(5) The intrinsic parameters of our desktop camera and the ODIS camera have been determined using 
our code. Our calibrated results have also been cross-validated using Microsoft code; 

(6) A new radial distortion model is proposed so that the radial undistortion can be performed 
analytically with no numerical iteration; 
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(7) Based on the results of this work, some new apphcation possibihties have been suggested for our 
mobile robots, such as ODIS. 

The rest of the report is arranged as follows. First, some notations and preliminaries are given, such 
as camera pinhole model, intrinsic parameters, and extrinsic parameters. Then, the calibration method 
proposed in [1, 2] is re-derived with a correction to a minor technical error in [1, 2]. Using this method, 
calibration results of 3 different cameras are presented. Finally, several issues are proposed for future 
investigations and possible applications of camera calibration are discussed. 

2 Camera Projection Model 

To use the information provided by a computer vision system, it is necessary to understand the geometric 
aspects of the imaging process, where the projection from 3-D world reference frame to image plane 
(2-D) causes direct depth information to be lost so that each point on the image plane corresponds to 
a ray in the 3-D space [7]. The most common geometric model of an intensity imaging camera is the 
perspective or pinhole model (Figure 1). The model consists of the image plane and a 3-D point Oc, 
called the center or focus of projection. The distance between the image plane and Oc is called the focal 
length and the line through Oc and perpendicular to the image plane is the optical axis. The intersection 
between the image plane and the optical axis is called the principle point or the image center. As shown 
in Figure 1, the image of P"^ is the point at which the straight line through Oc and P'^ intersects the 
image plane. The basic perspective projection [8] in the camera frame is 





_ f 




yC 







(1) 



where P"^ = [X*^, Z'^]-^ is a 3-D point in the camera frame and p = \x'^ ^y'^Y is its projection in the 
camera frame. In the camera frame, the third component of an image point is always equal to the focal 
length /. For this reason, we can write p = [x'^, y'^]^ instead of p = [x*^, y'^, f]^. 




Figure 1: The perspective camera model 



3 Aspects that Real Cameras Deviate from Pinhole Model 

A real camera deviates from the pinhole model in several aspects. The most significant effect is lens 
distortion. Because various constraints in the lens manufacturing process, straight lines in the world 
imaged through real lenses generally become somewhat curved in the image plane. However, this 
distortion is almost always radially symmetric and is referred to as the radial distortion. The radial 
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distortion that causes the image to bulge toward the center is called the barrel distortion, and distortion 

that causes the image to shrink toward the center is called the pincushion distortion (See Figure 2). 
The center of the distortions is usually consistent with the image center. 





1 1 
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Barrel Distorbon Hticushioa Distortion 



Figure 2: The barrel distortion and the pincushion distortion 



The second deviation is the flatness of the imaging media. However, digital cameras, which have 
precisely flat and rectilinear imaging arrays, are not generally susceptible to this kind of distortion. 

Another deviation is that the imaged rays do not necessarily intersect at a point, which means there 
is not a mathematically precise principle point as illustrated in Figure 3. This effect is most noticeable 
in extreme wide-angle lenses. But the locus of convergence is almost small enough to be treated as a 
point especially when the objects being imaged are large with respect to the locus of convergence. 



Incoming Rays 




Incoming Rays 




Pinhole Model 



Figure 3: Lens camera deviates from the pinhole model in locus of convergence 



4 Camera Parameters 



Definition: Camera Parameters 

Camera parameters are the parameters linking the coordinates of points in 3-D space with the coordi- 
nates of their corresponding image points. In particular, the extrinsic parameters are the parameters 
that define the location and orientation of the camera reference frame with respect to the world refer- 
ence frame and the intrinsic parameters are the parameters necessary to link the pixel coordinates of 
an image point with the corresponding coordinates in the camera reference frame. 



4.1 Extrinsic Pcirameters 

The extrinsic parameters are defined as any set of geometric parameters that uniquely define the trans- 
formation between the world reference frame and the camera frame. A typical choice for describing 
the transformation is to use a 3 x 1 vector t and a 3 x 3 orthogonal rotation matrix R such that 
P"^ = RP"' -I- t. According to Euler's rotation theorem, an arbitrary rotation can be described by 
only three parameters. As a result, the rotation matrix R has 3 dcgrcc-of- freedom and the extrinsic 
parameters totally have 6 degree of freedom. Given a rotation matrix R in Equation (2), one method 
to get the 3 parameters that uniquely describe this matrix is to extract ZYZ Euler angles [9] , denoted 
by (a, b, c), such that 
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r23 
r33 



where 



R,(c) 



cos(c) 
sin(c) 




R = R^(a) Kyib) R,{c), 

Ryib) 



- sin(c) 
cos(c) 




cos (6) sin (6) 

1 
— sin(6) cos(6) 



When sin{b) ^ 0, the solutions for (a, b, c) are 

b = arctan2 (y^ 



31 + ^32' ^33), 

a = arctan2 (r23/sin(6),ri3/sin(6)), 
c = arctan2 (r32/sin(6), — r3i/sin(6)). 



(2) 
(3) 

(4) 



(5) 



4.2 Intrinsic Parameters 

The intrinsic parameters cirG clS follows: 

• The focal length: / 

• The parameters defining the transformation between the camera frame and the image plane 
Neglecting any geometric distortion and with the assumption that the CCD array is made of 
rectangular grid of photosensitive elements, we have: 

y" = -{v- vo) Sy 

with {uo, vq) the coordinates in pixel of the image center and Sx, Sy the effective sizes of the pixel 
in the horizontal and vertical direction respectively. Let fx = f /Sx, fy = f /Sy, the current set of 
intrinsic parameters are uo,vo,fx, and fy. 

• The parameter describing the skewncss of the two image axes: J = fy tan 9 
The skewness of two image axes is illustrated in Figure 4. 




Figure 4: Skewness of two image axes 



• The parameters characterizing the radial distortion: ki and k2 
The radial distortion is governed by the Equation [10] 

F{r) = r f{r) = r (1 + fcir^ + A;2r^ + fcgr^ + •••)• (7) 

Two coefficients for distortion are usually enough. The relationship between the distorted and the 
undistorted image points can be approximated using 

Ud = u + {u- Uo) [ki{x^ + y2) + k2ix^ + y^)^] 

Vd = v + {v- vo) [ki{x'^ + y2) + k2{x^ + y^)^] ^ ^ 

where {ud,Vd) are the real observed distorted image points and {u,v) the ideal projected undis- 
torted image points. So, till now, the set of intrinsic parameters are uq, vo, fx, fy, 7, ^i, and k2- 
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4.3 Projection Matrix 



With homogeneous transform and the camera parameters, we can have a 3 x 4 matrix M, called the 
projection matrix, that directly links a point in the 3-D world reference frame to its projection in the 
image plane. That is: 



M 



1 



A[R tl 



1 



a 7 uo 
1 



\R tl 



1 



(9) 



where A is an arbitrary scaling factor and the matrix A fully depends on the intrinsic parameters. The 

calibration method used in this work is to first estimate the projection matrix and then use the absolute 
conic to estimate the intrinsic parameters [1, 2]. From Equation (1) and (6), we have 



u = 



(10) 



Prom Equation (9), we have 



u 



with scaling factor \ = Z'^. Prom the above two equations, we get a = —f/sx = 
manner, (3 = —fy. 



(11) 

—fx- In a same 



5 Extraction of Feature Locations 
5.1 Calibration Object 

The calibration method illustrated here is a self-calibration method, which uses a planar calibration 
object shown in Pigure 5, where 64 squares are separated evenly and the side of each square is 1.3 cm. 



Figure 5: Current calibration object 

The procedures to extract the feature locations of the above calibration object are illustrated in 
Table 2. The input image is an intensity image. After thresholding it with a certain value (which is 
150 in our case), we can get a binary image. The binary image then goes through some Connected 
Component Labeling algorithm [11] [12] that outputs a region map where each class of the connected 
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pixels is given a unique label. For every class in the region map, we need to know whether or not it can 
be a square box. In our approach, this is done by first detecting the edge of each class and then finding 
the number of partitions of the edge points. If the number of partitions is not equal to 4, which means 
it is not a 4-sided polygon, we will bypass this class. Otherwise, we will fit a line using all the edge 
points that lie between each two adjacent partition points and thus get 4 line fits. The final output of 
this class is the intersections of these 4 lines that approximate the 4 corners of each box. After running 
through all the classes in the region map, if the number of detected boxes equals to the actual number 
of boxes in calibration object, we will record all the detected corners and arrange them in the same 
order as for the 3-D points in space (for a given calibration object, assume Z"^ = 0, we know the exact 
coordinates of the feature points in world reference frame and we need to arrange these feature points 
in certain order so that after detecting feature points in the observed images, we can have an algorithm 
to seek the map from a point in the world frame to its corresponding projection in the image plane). 
After detecting several this kind of images, we are fully prepared to do calibration calculation. 

Table 2: Procedures to Extract Feature Locations for One Input Image 

Threshold input intensity image (PGM) to make it binary (the threshold is 150) 
Find connected components using 8-connectivity method 

Loop for every class in the region map 

Select the class whose area is < 3000 and > 20 
Binary edge detection of this class 
Find partitions of the edge points 
If # of partitions = 4 

Line fit between each two adjacent partition points 
Output 4 line intersections 
End if 
End loop 

If the total # of intersections = 4 x number jo j J30xes An jcalibraticm job ject 

Arrange intersections in the same order as points in the world reference frame 
End if 



5.2 Bineiry Image Edge Detection (Boundciry Finding) 

A boundary point of an object in a binary image is a point whose 4-neighborhood or 8-neighborhood 
intersects the object and its complement. Boundaries for binary images are classified by their connec- 
tivity and by whether they lie within the object or its complement. The four classifications axe: interior 
or exterior 8-boTindaries and interior or exterior 4-boundaries [13]. In our approach, we use interior 
8-boundary operator, shown in Figure 6, which is denoted as: 

b = (1- (aoiV)) a, (12) 

where 

(1) a is the input binary image 

(2) b is the output boundary binary image 

(3) N is the 4-neighborhood: N{p{u, v)) = {y : y = {u ± j,v) or y = {u, v £ {0, 1}} 

(4) For each pixel p{u, v) in a, p{u, v) o N = minimum pixel value around p{u, v) in the sense of 
4-neighborhood 
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iBoxes + Interior Boundary; iZZ% 23) -> [ODOOOO] ^ 



Figure 6: Objects and their interior 8-boundary 
5.3 Partitions of Edge Points 

Given a set of points that characterize the boundary of some object, a common question is what shape 
this object is, when we try to use polygons, especially the convex polygons, to denote objects in the 
real world. The set of points can be the output of some range finding sensors, such as laser and sonar. 
Or, it can come from images captured by a camera and is preprocessed by some edge detector, which 
is just the case we are discussing. In our problem, we know beforehand that the region of interest is a 
square and we can use the scan line approximation method [14, 15] to find the number of partitions. 
The scan line approximation algorithm is described in Table 3. Figure 7 is an illustration. 

Table 3: Scan Line Approximation Algorithm [14, 15] 



Problem Definition 
Assumption: Object is described using a convex polygon 
Given: A set of data points that have already been sorted in certain order 
Find: Partition points 

Algoritlim 

Scan_Line_Approximation (startJndex, end_index, data_points) 
Draw a line connecting start point and ending point 

Calculate the maximum distance each point that lies between start_index and end_index to this line 
If the maximum distance is greater than a predefined threshold 

Record the index corresponding to the point that gives the maximum distance 

Set end_index = the index of that point that gives the maximum distance 
Scan_Line_Approximation (startJndex, end judex, data_points) 

Set startJndex = the index of that point that gives the maximum distance 
Scan_Line_Approximation (startJndex, endJndex, data_points) 
End if 



In Table 3, the algorithm is described/implemented in a recursive way. Applying this algorithm, an 
important issue is how to decide the threshold. Unfortunately, this threshold is application-related. In 
our implementation, we choose 5 — 10 pixels. The smaller the boxes or the farther that the camera is 
way from the calibration object, the smaller the threshold should be. Figure 8 shows the fitting results 
using the partitions found by scan line approximation algorithm, where all the input data are the edge 
points of some classes in the region map. Another thing that we need to pay attention to is about 
how to choose the initial starting and ending points. It is obvious that they cannot be on the same 
side. Otherwise, due to noise in the data, the point whose distance to the line connecting starting and 
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Figure 7: Illustration of scan line approximation algorithm (a) 3 sides (b) 4 sides 

ending points is maximal might not be around corners. That is why we always start around corners as 
in Figure 7. This problem can be solved simply by first finding the two adjacent points whose maximal 
distance that all other points to this line is the biggest. 




Figure 8: Fitting results using partitions found by scan line approximation algorithm 

Figure 9shows an example of the processed images at all steps, where the input images are captured 
by a desktop camera. Notice that in Figure ??, in the image titled with "Binary Image + Partition 
Points", the triangular in the upper right corner does not show in the next step. The reason why this 
happens is that after the process of finding partition points, the number of partition points does not 
equal to 4 and we thus bypass this class. 

6 Calibration Method 

In this section, the calibration method in [1, 2] is described in detail. Using the calibration object shown 
in Figure 5, this algorithm is a self-calibration method. It only requires the camera to observe a planar 
pattern at a few different orientations. Either the camera or the calibration object can be moved by 
hand and the motion need not be known. The reason why this is feasible is that one image observed 
by a camera can provide 2 constraints about this camera's intrinsic parameters that are regarded to be 
unchanged here. With 3 images observed by the same camera, 6 constraints are established and we are 
able to recover the 5 intrinsic parameters. Once the intrinsic parameters are known, we can estimate the 
extrinsic parameters, the distortion coefficients {ki,k2), and put every initial guess of these parameters 
into some nonlinear optimization routine to get the final estimations. Another aspect that makes [1, 2] 
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Figure 9: Feature points extraction for desktop images (1) 



appealing is that the author provides cahbration results and an executable file on the web page [16] 
along with the sample images. The procedures to do calibration are illustrated in Table 4. 

The idea to assume the calibration object is always at Z"^ = even after some unknown movement 
maybe bewildering (We are talking about the case when we keep the camera static and move the 
calibration object). The common sense about the world reference frame is that it is unique. So, how 
can we assume the calibration object is still at Z"^ = after some rotation and translation? The 
answer to this question is: as mentioned before, only the relative position and orientation between 
the calibration object and the camera is of concern. Each image can provide 2 constraints that are 
independent to all others. Thinking in the other way, it is the same when we keep the calibration object 
static and move the camera. The basic calibration equations are given as follows. 



6.1 Homography Between the Model Plane and Its Image 

Without loss of generality, we assume the calibration object is Z^ = in the world reference frame. 
Let's denote the i^^ column of the rotation matrix R by ri, we have: 



M 



1 



A[ri r2 ra t] 



1 



A[ri r2 t 



1 



(13) 



Therefore a model points in 3-D space is related to its image by a homography H 





u 






A 


V 


= H 






1 




1 



(14) 



where 



H = A[ri ra t]. 

In this way, the 3x3 matrix H is defined up to a scaling factor. 

Given an image of the calibration object, the homography H can be estimated by maximum like- 
lihood criterion. Let Mi and rrii be the model and its image point respectively. Let's assume mi is 
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Table 4: Camera Calibration Procedures 



Linear Parameter Estimation 
Estimate Homographies (Section 6.1) 
Let be the number of images that we want to observe 
Loop for i from 1 to A'^ 

Assume the calibration object is at = 

Establishes the i*'* homography between the calibration object and its image 
Change the orientation of either calibration object or camera 
End Loop 

Estimate Intrinsic Parameters (Section 6.3) 

For each homography we have 2 constraints concerning the 5 intrinsic parameters 
Now we have 2N constraints and we can solve the 5 intrinsic parameters using SVD 

Estimate Extrinsic Parameters (Section 6.4) 

Using the estimated intrinsic parameters and homographies, we can estimate the extrinsic parameters 
Estimate Distortion Coefficients (Section 6.5) 

Using the estimated intrinsic and extrinsic parameters, we can get the ideal projected image points 
Along with the real observed image points, we can estimate the two distortion coefficients {ki,k2) 

Nonlinear Optimization 

(Section 6.6) 

Take all parameters estimated above as an initial guess 

Use some nonlinear optimization routine, we can get the final estimated values 



corrupted by Gaussian noise with mean and covariance matrix Aj„.. 
estimation of H is obtained by minimizing 



mi) 



Then the maximum likelihood 



(15) 



where h,- is the i*'* row of H and 



m,; = 



(16) 



In practice, we simply assume A^. = cr^ I for all i. This is reasonable if the points are extracted 
independently with the same procedure. For each pair of Mj and m^, we have 



u = 



Let X = [hj", h^^, hg'], then 



h^ilX^ + h2.2Y'" + h33 
h2lX'" + h22Y'^ + h23 

hsiX^ + h32Y'^ + hss' 



(17) 











1 
X"^ 







-uX"" -uY"" -U 
-vX"" -vY"^ -V 



(18) 
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When we are given n points, we have n above equations and we can write them in matrix form as 
Lx = where 



Xf 






yf 1 

xf yf 1 

y2"' 1 

1 



-uiXf 

-vixr 

-U2Xf 
-V2Xf 



-V2Y^ 



-Ul 
-Vl 
-U2 
-V2 



X- y» 1 

Xf 





1 



^ n 



-V X^ 

-V X^ 



The matrix L is a 2n x 9 matrix and the solution is well known to be the right singular vector of L 
associated with the smallest singular value. 

6.2 Constraints on the Intrinsic Parameters 
Given the estimated homography H = [hi,h2,h3], we have 

[hi,h2,h3] = A A [ri,r2,t], 
with A an arbitrary scalar. Using the knowledge that ri,r2 are orthogonal, we have 





(19) 



T 

ri r2 
T 

ri ri 



T 

r2 r2. 



Since Ari = A ^hi, Ar2 = A ^h2, 



hf A-^A-ih2 = 

hf A-^A-ihi = h^A-^A-ih2. 

Given a homography, these are the 2 constraints we obtained on the intrinsic parameters. 

6.3 Estimation of Intrinsic Pcirameters 

Let 



(20) 



(21) 



B = A-^A"^ = 



2_ 



Bu Bi2 5i3 

B21 B22 B23 

B31 -B32 -B33 

_ T 



■y(voj-uol3) 



/32 



7(t)o7-Mo/3) 
Q.2^2 



l_ -L 1 



(22) 



Note that B is symmetric. Define b = [511,^12,-622,-613,523,533]-^. Let the i^^ column vector of H 
be hi = [hii,hi2, ^is]^, we have 



hfBh, 



Denote 



= [hii hi2 his] 



Bii B12 5i3 

521 522 523 

531 532 533 
= hji{hiiBii + hi2Bi2 + hi^Bis 
+ hj3{hiiBi3 + hi2B23 + hisB^s). 











h,2 










+ hj2{hiiBi2 + hi2B22 + hi3B23) 



Vjj = [/iji^ji, hiihj2 + hi2hji, hi2hj2, hishji + huhjs, hizhj2 + hi2hj3, hishjs]'^ , 



(23) 
(24) 
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the two constraints in Equation (21) become 

■T 1 

• = 0. (25) 



(Vii - Vi2)^ 

If N images of the cahbration object are taken, by stacking N such equations, we have Vb = where 
V is a 2N X 6 matrix. 

When A'" > 3, we will have a unique solution defined up to a scaling factor. The solution is well 
known to be the right singular vector of V associated with the smallest singular value. The matrix B 
is estimated up to a scaling factor B = A A~^A~^. After estimation of B, the intrinsic parameters can 
be extracted from B by 

^^0 = (-B12-B13 — BiiB2^) / {B11B22 - -8^2)) 

A = -B33 — [i?i3 + Uo(-Bl2-Bl3 — i?ii-B23)]/-Bll, 



a = JX/Bn, (26) 



p = ^JXBu/{BuB22-Bl^), 

7 = -Bi2a^(3/X, 
uq = jvo/(3 - Bisa'^/X. 

The original equation to estimate uq in [1, 2] is uq = ^vo/a — Bi^a^/X. This must be an obvious 
mistake since when all the other 5 parameters are known, uq can be estimated directly from B13. The 
reason why using a wrong equation to estimate uq still achieves a good accuracy might due to the fact 
that a and (3 are the scaling factors in the two image axes and they are usually close to each other. 

6.4 Estimation of Extrinsic Pcirameters 

Once A is known, the extrinsic parameters can be estimated as: 

ri = AA-^hi 

r2 = AA-^ha (27) 
ra = ri X r2 

13 



t = AA-^ha 



where A = l/||A^^hi||2 = l/||A^^h2||2- Of course, due to the noise in the data, the computed 
matrix R = [ri r2 r^] does not satisfy the properties of a rotation matrix RR-^ = I. One way to 
estimate the best rotation matrix from a general 3x3 matrix R is: by the Matlab function svd with 
[U, S, V] = svd (R). The best rotation matrix will be UV^. 

6.5 Estimation of Distortion Coefficients 

Assume the center of distortion is the same as the principal point. Equation (8) describes the relationship 
between the ideal projected undistorted image points (n, v) and the real observed distorted image points 
{ud, Vd). Given n points in N images, we can stack all equations together to obtain totally 2Nn equations 
in matrix form as Dk = d, where k = [ki, /c2]"^. The linear least-square solutions for k is 

k = (D^D)-^D^d. (28) 

6.6 Nonlinecir Optimization: Complete Mctximum Likelihood Estimation 

Assume that the image points are corrupted independently by identically distributed noise, the maxi- 
mum likelihood estimation can be obtained by minimizing the following objective function 

N n 

J = Y1 Y.W'^ii - m(A, fci, fc2, Ri, ti, Mj)||2, (29) 
i=i j=i 
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where m(A, ki, A;2, Rj, tj, Mj) is the projection of point Mj in the i image using the estimated param- 
eters. This is a nonhnear optimization problem that can be solved by Matlab optimization function 
fminunc. 

In our implementation, one observation is that without an initial estimation of distortion coefficients 
{ki, ^2)) simply setting them to be 0, gives the same optimization results as the case with a good initial 
guess. Clearly, a "good" initial guess of the distortion coefficients is not practically required. 



7 Calibration of Different Cameras 

In this section, some calibration results are presented using the images provided in [16] . Images captured 
by a desktop camera and the ODIS camera are also used. For each camera, 5 images are captured and 
the feature locations are extracted for each image. Here, we always use 5 images for calibration. Using 
different number of images is also feasible. In practice, we found that 5 images are sufficient for camera 
calibration. 

7.1 Code Validation Using Images in [1, 2] 

In [1, 2], the calibration images are posted on web page [16]. We use the reported results to validate 
our implementation code with the same calibration images. 

7.1.1 Plot of the Observed and the Projected Image Points - Microsoft Images 

Figure 10 shows the observed and the projected image points using the images provided in [1, 2]. 




Figure 10: Plot of the observed and the projected image points - Microsoft images 
^Blue dots are the real observed image points 

^Red dots are the projected image points using the estimated camera parameters 



7.1.2 Comparison of Calibration Results - Microsoft Images 

The parameters before and after nonlinear optimization along with the Microsoft calibration results, 
obtained by executing its executable file posted on the web page [16], are shown in Table 5. Comparing 
these final calibration results, one can find that there are slight differences between some of the pa- 
rameters such as vq. However, as can be seen in section 7.1.3, the objective functions from these two 
calibration codes are very close. 
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Table 5: Comparison of Calibration Results - Microsoft Images 





(Jur Implementation 


Microsort 




Before Upti 


After Upti 


After Upti 


a 


871.4450 


832.5010 


832.5 


7 


U.2419 


0.2046 


0.2045 


1*0 


300.7676 


303.9584 


303.959 




871.1251 


832.5309 


832.53 




2:20.8684 


206.5879 


206.585 


ki 


0.1371 


-0.2286 


-0.2286 


k2 


-2.0101 


0.1903 


0.1903 



7.1.3 Objective Function - Microsoft Images 

Table 6 shows the comparison of the final values of objective function defined in Equation (29) after 
nonlinear optimization between Microsoft result and our implementation. The results show that they 
are very close. We can conclude that our code is correct for the Microsoft images. In what follows, we 
shall present two more groups of calibration results for our desktop camera and the GDIS camera in our 
center. As in the case of Microsoft images, we will compare the results similarly for further validation. 

Table 6: Objective Function - Microsoft Images 



Microsoft 


Our Code 


144.8799 


144.8802 



Remark 1 The options when using the Matlab function fminunc is not recorded. So, when using 
different options, slightly different results can be achieved. 

7.1.4 Nonlinear Optimization Iterations - Microsoft Images 

Table 7 shows the nonlinear optimization iterations, where the initial guess of all parameters are the 
estimations obtained in Section 6. The nonlinear optimization using Matlab function fminunc. From 
this table, we can see that after 52 iterations, the value of f{x), the objective function defined in 
Equation (29), drops to 144.88 from the initial value of 1055.89. 
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Tabic 7: Nonlinear Optimization Iterations - Microsoft Images 



Iteration 


Function 




Step-size 


Directional 




Count 


Derivative 


1 


37 


1055.89 


0.001 


-5.05e+009 


2 


78 


1032.26 


9.36421e-009 


-7.97e+004 


3 


120 


915.579 


1.85567e-007 


-1.13e+004 


4 


161 


863.19 


1.1597e-008 


-2.78e+005 


5 


202 


860.131 


1.77145e-008 


-2.45e+004 


6 


244 


836.386 


1.37495e-007 


-4.83e+003 


7 


285 


820.765 


1.58388e-008 


-1.13e+005 


8 


327 


816.13 


1.86391e-007 


-2e+003 


9 


368 


800.842 


3.79421e-008 


-4.12e+004 


10 


410 


788.888 


3.68321e-007 


-3.29e+003 


11 


452 


769.459 


1.47794C-006 


-1.38e+003 


12 


493 


738.541 


1.13935C-006 


3.2e+003 


13 


535 


692.716 


3.81991e-006 


606 


14 


576 


674.548 


1.0489e-006 


360 


15 


618 


631.838 


3.81559C-006 


-1.26e+003 


16 


659 


616.973 


1.91963C-006 


-1.79e+003 


17 


700 


604.857 


2.96745e-006 


-724 


18 


742 


593.179 


6.19011e-006 


-980 


19 


783 


573.66 


3.68278e-006 


-4.67e-|-003 


20 


824 


544.497 


4.69657e-006 


-2.28e+003 


21 


865 


537.636 


6.60037C-006 


-15.1 


22 


906 


530.468 


3.62979C-006 


-103 


23 


947 


525.032 


2.35736C-006 


-85.4 


24 


989 


523.091 


4.80959C-006 


22.9 


25 


1031 


509.698 


2.98894e-005 


-320 


26 


1072 


505.972 


2.3338e-006 


-2.54e-f003 


27 


1114 


499.005 


6.5114e-005 


-61.8 


28 


1155 


493.817 


4.78348e-006 


-587 


29 


1197 


468.823 


0.000433333 


-162 


30 


1238 


378.238 


0.000463618 


-924 


31 


1279 


280.68 


0.00024013 


-1.45e+003 


32 


1320 


230.465 


0.000148949 


-990 


33 


1361 


223.328 


0.00017819 


31.5 


34 


1402 


221.525 


0.0649328 


10.6 


35 


1444 


218.413 


0.179589 


-0.000716 


36 


1486 


191.642 


0.577805 


-0.521 


37 


1528 


176.735 


1.68909 


0.003 


38 


1570 


173.092 


1.62593 


8.26C-005 


39 


1612 


169.211 


1.38363 


-0.0033 


40 


1654 


157.984 


2.9687 


0.00291 


41 


1696 


146.471 


1.87697 


-0.0366 


42 


1737 


145.015 


0.894999 


0.00491 


43 


1778 


144.911 


0.869282 


-0.000345 


44 


1820 


144.893 


1.56986 


0.000205 


45 


1862 


144.882 


1.50158 


5.98C-005 


46 


1903 


144.88 


1.31445 


-0.00126 


47 


1946 


144.88 


0.0356798 


0.000261 


48 


1989 


144.88 


1.06404 


-9.25e-006 


49 


2030 


144.88 


0.771663 


1.23C-005 


50 


2071 


144.88 


0.123009 


-0.000482 


51 


2109 


144.88 


0.0615043 


-0.000114 


52 


2147 


144.88 


-0.0307522 


-1.18e-005 
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7.2 Calibration of a Desktop Camera 



This section shows the cahbration results for a desktop camera. 



7.2.1 Extracted Corners in the Observed Images - The Desktop Camera Case 

Figure 11 shows the extracted feature points in the observed images captured by the desktop camera. 
The extracted corners are marked by a cross and the dot in the center of each box is just to test if 
the detected boxes are in the same order as the 3-D points in the world reference frame. Due to the 
low accuracy of this camera, the extracted feature points deviate a lot from their "true positions", as 
"sensed" or "perceived" by our human observers. 
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Figure 11: Extracted corners in the observed images captured by the desktop camera 



7.2.2 Plot of the Observed and the Projected Image Points - The Desktop Camera Case 

Figure 12 shows the observed and the projected image points captured by the desktop camera. For 
descriptions, please refer to Figure 10. 
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Figure 12: Plot of the observed and the projected image points - the desktop camera case 



7.2.3 Comparison of Calibration Results - The Desktop Camera Case 

Table 8 show the calibration results of our implementation and Microsoft executable file. 

7.2.4 Objective Function - The Desktop Camera Case 

From Table 9, we can see that the final calibration results by our implementation and by the Microsoft 
group are almost identical. 
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Table 8: Comparison of Calibration Results - The Desktop Camera Case 





Our Implementation 


Microsoft 




Before Opti 


After Opti 


After Opti 


a 


350.066701 


277.1457 


277.145 


1 


1.693062 


-0.5730 


-0.573223 




200.051398 


153.9923 


153.989 




342.500985 


270.5592 


270.558 




100.396596 


119.8090 


119.812 


h 


0.096819 


-0.3435 


-0.343527 


k2 


-0.722239 


0.1232 


0.123163 



Table 9: Objective Function - The Desktop Camera Case 



Microsoft 


Our Code 


778.9763 


778.9768 



7.2.5 Nonlinear Optimization Iterations - The Desktop Camera Case 

For data format and descriptions, please refer to Table 7. 
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Table 10: Nonlinear Optimization Iterations - The Desktop Camera Case 
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7.3 Calibration of the ODIS Camera 

Now, we want to calibrate the ODIS camera, a 1/4 inch color board camera with 2.9 mm Lens [17]. 

7.3.1 Extracted Corners in the Observed Images - The ODIS Camera Case 

Figure 13 shows the extracted feature points in the observed images captured by the ODIS camera. 




MM M m M m.m 

jHJiiipdBintt 
Mmm0 mm** 

J!L« J« K EJ/ 




M M * * m m m m. 

mm m urn*** 

H a art at jM 

::owJ^*•■aw«jlT^^^f■a>^■;^^-^Jiftf^a(l 



Figure 13: Extracted corners in observed images captured by the ODIS camera 

7.3.2 Plot of the Observed and the Projected Image Points - The ODIS Camera Case 

Figure 14 shows the observed and the projected image points captured by the ODIS camera. For 
descriptions, please refer to Figure 10. 
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Figure 14: Plot of the observed and the projected image points - the ODIS camera case 



7.3.3 Comparison of Calibration Results - The ODIS Camera Case 

(See Table 11) 

7.3.4 Objective Function - The ODIS Camera Case 

(See Table 12) 

7.3.5 Nonlinear Optimization Iterations - The ODIS Camera Case 

(See Table 13) 
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Table 11: Comparison of Calibration Results - The ODIS Camera Case 





Our Implementation 


Microsoft 




Before Opti 


After Opti 


After Opti 


a 


320.249458 


260.7636 


260.764 


7 


10.454189 


-0.2739 


-0.273923 


Mo 


164.735845 


140.0564 


140.056 




306.001053 


255.1465 


255.147 


Vo 


85.252209 


113.1723 


113.173 




0.071494 


-0.3554 


-0.355429 


k2 


-0.3428(i(j 


0.1()33 


0.163272 



Table 12: Objective Function - The ODIS Camera Case 



Microsoft 


Our Code 


840.2665 


840.2650 
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Table 13: Nonlinear Optimization Iterations - The ODIS Camera Case 
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8 Conclusions and Further Investigations 



In this report, we have documented a complete implementation of the flexible camera calibration method. 
Using 3 cameras, we have cross- validated that our camera calibration code is as good as the Microsoft 
code posted on the web page [16]. We have re-derived all the equations in [1, 2] and corrected a technical 
error found in [1, 2]. Compared to the work in [1, 2], where the feature location algorithms were not 
discussed, we have built a complete code for camera calibration starting from the raw image acquisition 
step. A new method to effectively find the feature locations of the calibration object has been used in 
the code. More specifically, the scan line approximation algorithm is proposed to accurately determine 
the partitions of a given set of points. 

Based on our own platform and code, we can try some new ideas. In what follows, we will describe 
some of the efforts for the improved performance in camera calibration. 

8.1 Error Reduction in Extraction of Feature Locations 

A big error source for calibration is the error in extraction of the feature locations and this is noticeable 

in Figure 11. There arc a number of available feature extraction procedures developed [18, 19] and 
the procedure in Section 5 is our local line fitting method. Since the line fitting method is based on 
a small number of pixels, each noisy data point will strongly affect the fitting result. To improve the 
accuracy, an alternate method for line fitting is to use the Hough Transform, which is a standard tool 
in the domain of artificial vision for the recognition of straight lines, circles [20] and ellipses [20] and is 
particularly robust to missing and contaminated data [21, 20]. 

Using the calibration object in Figure 5, one shortcoming is that when choosing different threshold 
values to convert the observed intensity image to the binary image, the threshold may cause the squares 
to shrink or enlarge, which will in turn affect the feature localization. In [6, 18], the checkerboard 
pattern is proposed as the calibration object. The advantage of the checkerboard pattern, as shown in 
Figure 15, is that its corners are localizable independent of the linearity of the image response. That 
is, applying a nonlinear monotonic function to the intensity values of the checkerboard image, such as 
gamma correction, does not affect the corner localizations. Using the checkerboard calibration pattern, 
the feature locations can be extracted by first convolving the mask image (obtained by thresholding the 
input PGM format image) with the window shown in Table 14 [18]. Since this filter itself resembles 
a checkerboard pattern, it gives a strong response, positive or negative, depending on which type of 
corner when centered over a checkerboard corner. The filter output of the mask image produces an 
image where the checkerboard corners appear as white or black dots (Sec Figure 16 and 17). In [18], 
localizing a particular checkerboard corner after the filter convolution is said to be easily accomplished 
by locating the point of maximum positive or negative filter response. Though not so straightforward, 
the checkerboard corners can be extracted by adding two extra steps. In the image of filter output, 
the values of most pixels are except in the small region around corners or edges, where pixels can 
be mostly positive or negative depending on what kind of corners. So, the first step added is to get 
the region map and in each region locate the pixel whose filter response is the strongest. Doing this 
way, the output pixels may include some noisy pixels around the edges and the second step is to pick 
up a certain number of pixels with the strongest filter response from all regions based on the known 
knowledge of how many corners in the given pattern. The filter output and extracted corners are shown 
in Figures 16 and 17. As can be seen from Figure 17, the current simple method is not robust enough 
and the feature locations can not always be extracted accurately. Some post-processing algorithms are 
thus needed. Other algorithms exist such as 1) using Canny edge detector to find the edge pixels, 2) 
globally fit curves to the edge pixels and 3) find their intersections [19] (also using the checkerboard 
pattern). However, for both calibration objects, in order to achieve a sub-pixel accuracy, a series of 
post-processing algorithms are necessary. 



22 



Figure 15: The checkerboard calibration object 



Table 14: Convolution Window for the Checkerboard Calibration Object 
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Figure 16: The result of convolving mask image with filter in Table 14 (1) 




Figure 17: The result of convolving mask image with filter in Table 14 (2) 
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8.2 Speed Up the Nonlinear Optimization 



Currently, we put totally 7 + 6A^ parameters into the nonlinear optimization routine to get the final 
estimations. In Matlab, this is absolutely time consuming. Can we make nonlinear optimization faster? 

Apply 5 images to get the initial guess of camera parameters. Then use 3, 4, and 5 image data for 
nonlinear optimization respectively. The 7 intrinsic parameters at each iteration are plotted in Figure 
18, where blue color is for the nonlinear optimization using 3 images, similarly red for 4 images, and 
green for 5 images. It is clear that the final estimation are very consistent with each other and using 
only 3 images for nonlinear optimization is absolutely time saving. 




Figure 18: Intrinsic parameters through nonlinear optimization iterations 
Blue color is for 3 images, red for 4 images, and green for 5 images 



8.2.1 Incremental Nonlinear Optimization 

To speed up the nonlinear optimization, one idea we came out with is to apply the first 3 images to 
nonlinear optimization routine, and then use the output 7 intrinsic parameters (treat these 7 parameters 
as unchanged) to optimize the extrinsic parameters for the 4*^^ and 5*'* image. At last, put all the 
estimated parameters again into the nonlinear optimization routine to get the final estimates. Doing 
this way, we call it incremental nonlinear optimization. Disappointedly, compared with the one-step 
optimization method, our initial attempt shows that the incremental method can only reach the same 
performance in terms of the accuracy and the consuming time. 

8.2.2 Constrained Nonlinear Optimization 

Another idea that might speed up the nonlinear optimization is to add some constrains, such as: 
• Mo is close to half of image's length 
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• vo is close to half of image's width 

• angle between two image axes is close to 90 degree 

• a, f3 arc greater than 

We shall investigate this idea later. 

8.3 Radial Undistortion 

According to the radial distortion model in Equation (8) , radial distortion can be understood to perform 
in one of the following two ways: 

• Perform distortion in the camera frame, then transform to the image plane 



X 




x' 














y 




y' 




Vd 



• Transform to the image plane, then perform distortion in the image plane 



X 




u 




Ud 


y 




V 




Vd 



To explain, let /(r) = 1 + kir^ + k2r^, = + y^, we know that 

Ud = {u- uo) f{r) + no 

= a x/(r) + 7 y/(r) +no 
= ax' + ^y' + uo, 
{v - -^o) f{r) + Vo 

P y' + vo, 



Vd 



where 



x' = X f{r) 
y' = y fir). 



(30) 
(31) 

(32) 



Radial undistortion is to extract {u,v) from {ud,Vd), which can be done by extracting {x,y) from 
{x\y'). The following derivation shows the problem when trying to extract (x,y) from (x',y') using 
two distortion coefficients ki and k2- From {ud,Vd), we can calculate {x',y') by 



(33) 
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Now the problem becomes to extract {x,y) from {x',y'). According to Equation (32), 

x' = x /(r) = a; [1 + ki{x^ + y^) + ^2(^2 + y^)^] 
y' = y f{r) = y [1 + ki{x'^ + y^) + k2{x'^ + y'^f]. 



(34) 



Let c = y'/x' = y/x, we have y = cx where c is a constant. Substituting y = cx into the above equation, 
gives 



x' = x[l + ki{x^ + c\^) + k2ix^ + c^x^f] 
= x + ki{l + c^)x^ + k2{l + c'^fx^. 



(35) 
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Let fi{x) = X + A;i(l + c^)x^ + k-zil + c^)^x^, then = -fi{x) so /i(.t) is an odd function. The 

analytical solution of Equation (35) is not a trivial task. It is said that this problem is still open (of 
course, we can use numerical method to solve it). But if we set k2 = 0, the analytical solution is 
available and radial undistortion can be done easily. In [22], for the same practical purpose, only one 
distortion coefficient ki is used to approximate radial distortion. 

Here, let's consider the distortion model again. There are two questions we are of particular interests. 
The first question is: can we choose 

radial distortion model case2 : /(r) = 1 + kir^ (36) 

instead of 

radial distortion model casei : /(r) = 1 + /cir^ + k2r'^ ? (37) 

What can be the performance degradation? The second question is: can we model radial distortion in 
other way so that we can achieve good and reasonable accuracy along with easy analytical undistortion? 
To answer the first question, we can re-optimize the parameters with only one distortion coefficient ki. 
Recall that the initial guess for radial distortion is done after having estimated all other parameters and 
just before nonlinear optimization. So we can reuse the estimated parameters and choose the initial 
guess for ki to be 0. For the second question, let's take 

radial distortion model cases • f{f) = 1 + kir + k^r^, (38) 

which is a function only related to radius r. The motivation of choosing this radial distortion model is 
that the resultant approximation of x' is also an odd function of x, as can be seen next. 
When F{r) = rf{r) = r(l + kir + k2r^), we have 



x' = X f{r) =x [(1 + kir + k2r^)] 
y' = y f{r) = y [(1 + hr + k2r'^)]. 



(39) 



Again let c = y' /x' = y/x, we have y = cx where c is a constant. Substituting y = cx into the above 
equation, gives 



x' = X 



X 



1 + ki Vx^ + C2x2 + k2ix'^ + C^X^] 



1 + ki\/l + c2 sign(x)x + ^2(1 + c^)x" 



= X + kWl + c2 sign(x) x^ + k2{l + c^) x^ . (40) 

Let /3(x) = X + fci Vl + sign(x) x^ + ^2(1 + c^) x^, fz{x) is also an odd function. 

Now, we want to compare the three radial distortion models based on the final value of objective 
function after nonlinear optimization by Matlab function fminunc. Using Microsoft images, desktop 
images, and ODIS images, the objective function, 5 intrinsic parameters, and distortion coefficients are 
shown in Table 15. The results in Table 15 show that the objective function of Models is always greater 
than that of Modeli, but smaller than that of Model2. This is consistent with our anticipation. 

The benefits of using radial distortion models are: 

(1) Low order fitting, better for fixed-point implementation 

(2) Explicit inverse function with no numerical iterations 

(3) Better accuracy than radial distortion model2 

To extract x from x' in Equation (40), let's first assume x > 0. Then there is an unique solution 
G IZ, denoted by x+ , that satisfies Equation (40) but may or may not coincide with the assumption. If 
x+ > 0, X = x_|_, otherwise, x = x_, which is the solution (g TZ) for the case when x < 0. 
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Table 15: Comparison of Three Distortion Models 





Microsoft Images 


Desktop Images 


ODIS Images 


Model 


#1 


#2 


#3 


#1 


#2 


#3 


#1 


#2 


#3 


J 


144.88 


148.279 


145.659 


778.9768 


904.68 


803.307 


840.2650 


933.098 


851.262 


a 


832.5010 


830.7340 


833.6623 


277.1457 


275.5959 


282.5664 


260.7636 


258.3206 


266.0861 


7 


0.2046 


0.2167 


0.2074 


-0.5730 


-0.6665 


-0.6201 


-0.2739 


-0.5166 


-0.3677 


Uo 


303.9584 


303.9583 


303.9771 


153.9923 


158.2014 


154.4891 


140.0564 


137.2155 


139.9177 


/3 


832.5309 


830.7898 


833.6982 


270.5592 


269.2307 


275.9040 


255.1465 


252.6869 


260.3145 




206.5879 


206.5692 


206.5520 


119.8090 


121.5254 


120.0952 


113.1723 


115.9295 


113.2417 


h 


-0.2286 


-0.1984 


-0.0215 


-0.3435 


-0.2765 


-0.1067 


-0.3554 


-0.2752 


-0.1192 


k2 


0.1903 





-0.1565 


0.1232 





-0.1577 


0.1633 





-0.1365 



8.4 Possible Applications of Camera Calibration in ODIS Missions 

8.4.1 Visual Servoing Continuity 

The wireless visual servoing in [23, 24] is performed by an uncalibrated camera. The experiment results 
do not match the theoretical simulation results satisfactorily because the parameters of the camera 
model are roughly set. With the camera calibration work presented here, we should be able to identify 
the parameters needed in simulation. 

8.4.2 Estimation of the Pan/Tilt Angles of the ODIS Camera 

The following steps can be applied to estimate the pan/tilt angles of the ODIS camera with respect to 
the ODIS body fixed coordinate system: 

• put the calibration object perpendicular to the ground, and let the bottom-left corner be [0, 0, 0] 

of the world coordinate system 

• put ODIS in front of the calibration object, and be sure to make its xb and ys axes parallel to 

those of world coordinate system 

• take a picture of calibration object using ODIS camera 

• estimate homography H 

• estimate 6 extrinsic parameters using intrinsic parameters we already know and the homography 
estimated in the last step 

• from the 3 angular vector, we can estimate the pan/tilt angles of ODIS camera 

8.4.3 Non-iterative Yellow Line Alignment with a Calibrated Camera 

Instead of our previous yellow line alignment method described in [23, 24], can we align to yellow line 
with a non-iterative way using a calibrated camera? The answer is yes. Let's begin with a case when 
only ODIS's yaw and x.y positions arc unknown while ODIS camera's pan/tilt angles are unchanged 
since calibration. The problem is described in detail in Table 16. 

Knowing that a change in ODIS's yaw angle only results in a change of angle c in the ZYZ Euler 
angles (a, b, c). So, when using ZYZ Euler angles to identify ODIS camera's orientation, we will assume 
the first two varialbes a, b are unchanged. In Figure 19, after some time of navigation, the robot thinks 
it is at Position 1. Then it sees the yellow line, whose locations in 3D world reference frame are known 
from map (denoted by and -Pg). After extracted the correspoinding points in image plane of the 
yellow line's two ending points, we can calculate the undistorted image points and thus recover the 
3D locations of the two ending points (denoted by and Pbb)i using ODIS camera's 5 intrinsic 
parameters and 2 radial distortion coefficients. Prom the difference between the yellow line's actual 
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Table 16: Task of Yellow Line Allignment Using Calibrated ODIS Camera 



Given: 3D locations of yellow line's two ending points 
ODIS camera's pan/tilt angles 
ODIS camera's intrinsic parameters 

ODIS camera's radial distortion model and distortion coefficeints 
The two projected points in image plane using ODIS camera 
Find: ODIS's actural yaw and x,y positions 



locations in map and the recovered locations, the deviation in the robot's x, y positions and yaw angle 
can be calculated. 
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Figure 19: The task of yellow line alignment 



Let {xi,yi) and {x2,y2) be the corresponding points in camera frame of yellow line's two ending 
points. Let R2 and t2 be the rotation matrix and translation vector at position 1 (where the vehicle 
thinks it is at), similarly Ri and ti at position 2 (the true position and orientation), we can write 
i?2 = ^R ■ Ri and t2 = ti + At, where AR and At are the deviation in orientation and translation. If 
the transform from the world reference frame to the camera frame is P'^ = R~^{P^ — t), first we will 
calculate P^^ and P^j^. Let PJ^ = [X]^^, FJ^, 0], we have 



pc 





— i?2 ^ 






Yaa - ^22 






— 123 



(41) 



Since 



yc 

yi 



T 



(42) 



we have two equations containing two variables and P^a calculated out. In a same way, we can 

get P^B- 

Knowing and Pg^, we have 



Xl 

yi 

1 



R2^AR{P^ - ti) = R^'iP^A - t2 



-1/ 



(43) 
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where A is a scaling factor. Prom Equation (43), we get 



i?2-1Ai?(Pl - ti) - + t2] = 0. 



(44) 



Similarly, we get 



R^'[AR{PS-ti 



)-P^S+t2]=0. 



(45) 



Using the above two equations, we have (-PXl ~ ^bb) ~ i^A ~ ^b)^ where /S.R is of the form 



So, AO is just the rotation angle from vector — > Pg to vector PJ^ — Pg^. Knowing AP, ti 
becomes ti = P^ - Ap-^P^^ + Ap-^P2t2- 

Remaining question: 

• Choose X from 6 possible values in radial undistortion for distortion Models 

• If GDIS camera's pan/tilt angles are changed due to user commanded tasks and are not available, 
what can we do for yellow line alignment? 

8.5 Tidy Up of Our Camera Calibration Code 

8.5.1 Singularity Reminder for Robustness in Camera Calibration 

When moving the camera or the calibration object to capture images, if the relative orientation between 
the calibration object and the camera docs not change, singularity will occur. We will add one function 
in the code that can pop up information telling the user which two images may cause the singularity 
problem. We will also design an automatic procedure to discard the image that cause the singularity 
problem and use the remaining images for calibration, if the number of remaining images is greater than 
2. To do so, a numerical indicator will be used to predict the degree of singularity. 

8.5.2 Two Complete Sets of Codes: Matlab and C/C++ 

Currently, all of the codes are in C/C++ except the nonlinear optimization routine. We should be able 
to find some open source code in C/C++ for nonlinear optimization so that wc can have a complete 
set of C/C++ code for camera calibration. Meanwhile, we would like to have a pure Matlab version by 
CMEX facility provided in Matlab. 

A Nonlinear Optimization using TOMLAB 

Using the nonlinear optimization problem defined in section 6.6 and with the data set for calibration of 
GDIS camera, wc can evaluate the performance of TOMLAB NLPLIB Toolbox with the results reported 
in the literature, such as Matlab Gptimization Toolbox. TGMLAB NLPLIB toolbox is a set of Matlab 
solvers, test problems, graphical and computational utilities for unconstrained and constrained opti- 
mization, quadratic programming, unconstrained and constrained nonlinear least squares, box-bounded 
global optimization, global mixed-integer nonlinear programming, and exponential sum model fitting 
[25, 26]. In TOMLAB, the routine ucSolve implements a prototype algorithm for unconstrained or 
constrained optimization with simple bounds on the variables, i.e. solves the problem 



AP 



cos(A6l) -sin(A6l) 
sin(A6') cos(A6l) 
1 



(46) 



miuj; f{x) 

S.t. XL < X < Xlf, 



(47) 
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Table 17: Algorithms Implemented in ucSolve 



Variable 


Algorithm 





Default algorithm (BFGS or Newton) 


1 


Newton with subspace minimization, using SVD 


2 


Safeguarded BFGS with standard inverse Hessian update 


3 


Safeguarded BFGS with Hessian update, and SVD or LU to solve 


4 


Safeguarded DFP with standard inverse Hessian update 


5 


Safeguarded DFP with Hessian update, and SVD or LU to solve 


6 


Fletcher-Reeves CG 


7 


Polak-Ribiere CG 


8 


Fletcher conjugate descent CG-method 



where x,xl, and xu G 7^" and /(x) G TZ. The algorithms implemented in ucSolve are listed in Table 
17. 

Performance comparison of these two nonlinear optimization toolboxes are based on the 
number of iterations and consuming time to reach the same or similar value for objective 
function. It should also be kept in mind that the comparison should go to the case when 
both these two toolboxes implement the same or similar algorithms. 

In Matlab, we choose function fminunc that uses BFGS Quasi-Newton method with a mixed 
quadratic and cubic line search procedure (because we do not provide gradient information and we 
thus use Medium-Scale Optimization algorithm) [27]. Correspondingly, in TOMLAB, we pick up al- 
gorithm #0 with Quadratic Interpolation or Cubic Interpolation LineSearchType. This is because 
ucSolve does not have a mixed quadratic and cubic line search procedure and the above two methods 
are what available. The command lines are: 

• Matlab fminvmc command line: 

options = optimset('Display','itcr', 'LargeScale', 'off', 'TolX', 10-^ 'TolFun', 10"^); 
X = fminunc (@ my fun, xO, options); 

Note: when 'LineSearchType" is not specified, the default value is "quadcubic" . 

• TOMLAB ucSolve command line: 

Prob = conAssign('myfun', [], [], [], [], [], 'calibration', xO); 
Prob.Solver.Alg = 0; 
Prob.PriLevOpt = 2; 

Result = ucSolve(Prob); 

Note: to specify the line search method for ucSolveo, parameter "LineAlg" is changed inside routine "LineSearch.m" . 

Table 18 and 19 show the optimization results for Matlab fminunc and TOMLAB ucSolve respec- 
tively, where each column represents Iteration Number, Objective Function at Each Iteration, Time at 
Each Iteration (obtained by Matlab command "rem(now,l)"), and Parameters Estimated respectively 
(in both tables, the number with a left arrow on the side is the total running time till current iteration). 
From these two tables, we observe: 

• Both functions converge 

• Both functions converge to close values for both objective function and estimated parameters 

• Matlab converge with less iterations and less time, slightly having advantage over TOMLAB 

For line search method, it is a common sense that quadratic interpolation that involves a data 
fit to the form ax^ + bx + c should take less time than cubic interpolation ax^ + bx^ + cx + d. 
In a same manner, a mixed quadcubic interpolation might lie somewhere between quadratic and 
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cubic interpolations. Thinking this way, Matlab function demonstrates faster convergent speed 
over TOMLAB. 

B Raw Data of the Extracted Feature Locations 

Please see corners.dat. 

C Matlab Source Code 

Please see codes.m. 
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Table 18: Matlab fminvmc with LineSearchType = QuadCubic 



Matlab fminunc 


# 


f{x) 


Time (10"^) 


Params 


1 


5872.79 


5.309561 


260.7636 


2 


5849.4 


5.310881 


-0.2739 


3 


5841.17 


5.312219 


140.0564 


4 


5734.46 


5.313533 


255.1465 


5 


5703.36 


5.314850 


113.1723 


6 


5524.51 


5.316187 


-0.3554 


7 


5505.25 


5.317472 


0.1633 


8 


5286.49 


5.318820 


-1.8822 


9 


5252.79 


5.320125 


2.8795 


10 


5147.55 


5.321412 


-1.9377 


11 


5123.77 


5.322762 


9.4005 


12 


5111.44 


5.324098 


-9.6023 


13 


5085.05 


5.325382 


-24.9324 


14 


5033.1 


5.326730 


-1.5881 


15 


5017.1 


5.328034 


2.8809 


16 


4875.16 


5.329351 


-1.6277 


17 


4638.53 


5.330698 


7.3082 


18 


4569.17 


5.332003 


-12.8265 


19 


4502.31 


5.333320 


-25.2393 


20 


4446.65 


5.334666 


-1.3550 


21 


4199.75 


5.335973 


2.8024 


22 


4140.98 


5.337258 


-1.3853 


23 


4052.51 


5.338576 


10.6469 


24 


3989.78 


5.339881 


-9.5356 


25 


3907.17 


5.341169 


-17.4115 


26 


3856.12 


5.342516 


-1.3737 


27 


3779.25 


5.343853 


2.8125 


28 


3645.96 


5.345139 


-1.3938 


29 


3619.6 


5.346456 


10.2720 


30 


3539.96 


5.347795 


-10.2214 


31 


3204.79 


5.349080 


-19.7591 


32 


3173.84 


5.350395 


-0.8848 


33 


3039.27 


5.351733 


2.6428 


34 


2869.31 


5.353052 


-0.8422 


35 


2805.64 


5.354367 


8.8168 


36 


2252.13 


5.355673 


-9.1196 


37 


1774.1 


5.356961 


-17.3228 


38 


1438.04 


5.358281 




39 


1242.37 


5.359631 




40 


1166.62 


5.360914 




41 


1124.36 


5.362229 




42 


1030.76 


5.363537 




43 


957.955 


5.364837 




44 


905.386 


5.366189 




45 


877.988 


5.367488 




46 


856.077 


5.368767 




47 


847.3 


5.370106 




48 


843.999 


5.371396 




49 


842.438 


5.372726 




50 


841.821 


5.374871 




51 


841.286 


5.376304 




52 


840.828 


5.377721 




53 


840.559 


5.379342 




54 


840.414 


5.380722 


^0.0071161 


55 


840.334 


5.382321 




56 


840.311 


5.383757 




57 


840.29 


5.385097 




58 


840.274 


5.386385 




59 


840.269 


5.387730 




60 


840.266 


5.389072 




61 


840.265 


5.390453 




62 


840.265 


5.391947 


^0.0082386 
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Tabic 19: TOMLAB ucSolve Algorithm, 



TOMLAB ucSolveo 


LineSearchType = Quadratic 


Interpolation 


LineSearchType = Cubic 


Interpolation 


Tr 


f(x) 


Time (10"^) 


Pcirams 


f(x) 


Time (10~^) 


Params 





5872.79 


7.905221 


260.5232 


5872.79 


5.607743 


260.6613 


1 


5850.93 


7.906629 


-0.2585 


5850.93 


5.614631 


-0.2650 


2 


5842.7 


7.908060 


139.4244 


5842.81 


5.623759 


139.4470 


3 


5736.13 


7.909437 


254.8903 


5750.96 


5.630199 


255.0380 


4 


5705.03 


7.910831 


113.0924 


5728.82 


5.637519 


113.0792 


5 


5589.15 


7.912188 


-0.3549 


5616.05 


5.642693 


-0.3553 


6 


5573.78 


7.913598 


0.1630 


5598.15 


5.650420 


0.1635 


7 


5370.21 


7.914950 


-1.8919 


5406.89 


5.656750 


-1.8910 


8 


5306.06 


7.916366 


2.8786 


5365.39 


5.663024 


2.8785 


9 


5231.59 


7.917751 


-1.9471 


5287.34 


5.669369 


-1.9462 


10 


5208.66 


7.919179 


9.3391 


5271.08 


5.678234 


9.3412 


11 


5195.82 


7.920595 


-9.6113 


5255.9 


5.685982 


-9.6124 


12 


5170.35 


7.921982 


-24.9336 


5237.29 


5.693702 


-24.9434 


13 


5165.73 


7.923401 


-1.5973 


5232 


5.702844 


-1.5973 


14 


5110.57 


7.924798 


2.8806 


5186.13 


5.712348 


2.8805 


15 


5050.42 


7.926171 


-1.6368 


5050.87 


5.721648 


-1.6368 


16 


4878.18 


7.927546 


7.2445 


4884.18 


5.728177 


7.2466 


17 


4796.39 


7.928953 


-12.8362 


4761.34 


5.737326 


-12.8376 


18 


4735.7 


7.930365 


-25.2287 


4705.51 


5.743843 


-25.2428 


19 


4688.86 


7.931775 


-1.3620 


4662.16 


5.751742 


-1.3618 


20 


4609.26 


7.933138 


2.8029 


4530.53 


5.760916 


2.8027 


21 


4505.35 


7.934551 


-1.3920 


4465.29 


5.767511 


-1.3918 


22 


4466.61 


7.935955 


10.6028 


4379.13 


5.775407 


10.6040 


23 


4373.21 


7.937358 


-9.5433 


4323.05 


5.783387 


-9.5435 


24 


4266.42 


7.938775 


-17.4184 


4191.96 


5.792649 


-17.4255 


25 


4220.26 


7.940171 


-1.3810 


4136.08 


5.803394 


-1.3805 


26 


4050.48 


7.941566 


2.8130 


4044.95 


5.808722 


2.8127 


27 


4023.58 


7.942962 


-1.4007 


3938.68 


5.815519 


-1.4002 


28 


3992.68 


7.944357 


10.2221 


3916.47 


5.824949 


10.2234 


29 


3767.4 


7.945721 


-10.2300 


3846.8 


5.831830 


-10.2301 


30 


3709.22 


7.947147 


-19.7626 


3576 


5.840273 


-19.7701 


31 


3614.32 


7.948520 


-0.8885 


3536.47 


5.849072 


-0.8886 


32 


3405.79 


7.949917 


2.6443 


3349.57 


5.856030 


2.6440 


33 


3191.31 


7.951232 


-0.8455 


3205.81 


5.860392 


-0.8456 


34 


3092.14 


7.952565 


8.7725 


3169.47 


5.864614 


8.7735 


35 


2669.06 


7.953911 


-9.1266 


2766.8 


5.867528 


-9.1268 


36 


1935.55 


7.955214 


-17.3202 


1872.21 


5.868979 


-17.3295 


37 


1856.89 


7.956483 




1534.06 


5.871888 




38 


1406.16 


7.957785 




1422.71 


5.873300 




39 


1215.79 


7.959123 




1247.39 


5.874730 




40 


1146.05 


7.960446 




1163.82 


5.877753 




41 


1062.19 


7.961723 




1105.05 


5.879195 




42 


990.388 


7.963026 




1035.38 


5.880571 




43 


944.654 


7.964327 




970.948 


5.881957 




44 


896.425 


7.965627 




916.323 


5.883337 




45 


868.715 


7.966990 




886.546 


5.884838 




46 


854.028 


7.969591 




864.105 


5.886252 




47 


846.043 


7.970898 




853.08 


5.889118 




48 


843.861 


7.972203 




848.483 


5.890505 




49 


842.991 


7.973501 




844.994 


5.891939 




50 


841.9 


7.974808 




843.241 


5.893339 




51 


841.566 


7.976111 




842.484 


5.894753 




52 


841.119 


7.977419 




841.523 


5.896165 




53 


840.647 


7.978721 




840.877 


5.897563 




54 


840.618 


7.980032 




840.686 


5.898996 




55 


840.558 


7.981342 




840.584 


5.900508 




56 


840.542 


7.982662 




840.515 


5.901898 




57 


840.542 


7.984747 




840.511 


5.903295 




58 


840.542 


7.985218 


<-0.0079997 


840.511 


5.914601 


<-0.0306858 
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